# Ryan Brutger and Alexandra Guisinger
# Replicatin code for Labor Market Volatility, Gender, and Trade Preferences
# Analysis run on macOS Mojave, MacBook Pro, 2.8 GHz Intel Core i7
# using R version 3.6.1 (2019-07-05) -- "Action of the Toes"

# ----------------------------------------------------------------------
# load all relevant themese and packages
# ----------------------------------------------------------------------

library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("RItools")
library("estimatr")
library("fastDummies")
library("ggpubr")
library("foreign")
library("mediation")
library("demoGraphic")

#Set your working directory
setwd() 

# Read data
data <- read_csv("Replication_JEPS_Volatility_Trade.csv")

# Calculate baseline employment prospects and trade support for men/women in control condition
data$enhanceEmp <- ifelse(data$employment==3, 1, 0) #indicator for belief that trade enhances employment opportunity
data$TRsupp<-ifelse(data$trade_support>3, 1, 0)#indicator for those who somewhat or strongly favor trade

# Calculate baseline support for trade and employment beliefs 
# for men and women in the control condition, and caonfidence intervals from model
trade.model <- lm(TRsupp ~ women, data=data[data$volatility==0, ])
summary(trade.model) #64% men support trade; .642-.125 = 51.7% women support trade (diff. p<0.01)
confint(trade.model) #95% CI: Men=(.616, .668); Women= (.481, .553)
.6417-.1606 #.481
.6417-.0888 #.553
# 
empl.model <- lm(enhanceEmp ~ women, data=data[data$volatility==0, ])
summary(empl.model) # 8.9% (p<0.01) more men believe trade enhances employment prospects
confint(empl.model) #95% CI: Men=(.262, .306); Women= (0.164, 0.226)
.284-.120 #0.164
.284-.058 #0.226

# Generate Figure 1 of the paper
base.names <- c("Men", "Women", "Men", "Women")
base.opin <- c(.642, .517, .284, .195)
color <- c("gray", "lightgray", "gray", "lightgray")

barplot(base.opin, main="", xlab="", ylim=c(0,.9), 
        ylab="Proportion of Respondents", names.arg= base.names, col=color, space=c(.2, .2, .8, .2))
#text(2.8, .85, labels = "Baseline Opinions in Control Condition", cex=1.5)
text(1.3, .7, labels = "Support Trade", cex=1)
text(4.3, .36, labels = "Trade Enhances \n Employment Propsects", cex=1) # export as 5.5 x 6.5 dimensions
lines(c(.7,.7), c(.616, .668), lwd = 4)
lines(c(1.9, 1.9), c(.481, .553), lwd = 4)
lines(c(3.72, 3.72), c(.262, .306), lwd = 4)
lines(c(4.92, 4.92), c(0.164, 0.226), lwd = 4)

# create data subsets data to only those who are working (working full time AND working part time) and not working 
data.working <- subset(data[data$workforce==1 | data$workforce==2 | data$workforce==3, ])
data.not.working <- subset(data[data$workforce!=1 & data$workforce!=2 & data$workforce!=3, ])

#Main Paper Table 1:
#Test H1a
vol.employment <- lm(employment ~ volatility, data = data.working)
#Test H1b
vol.lose.job <- lm(lose_job ~ volatility, data = data.working)
#Test H1c
vol.trade.benefits.self <- lm(trade_benefits.self ~ volatility, data = data.working)
vol.trade.benefits.region <- lm(trade_benefits.region ~ volatility, data = data.working)
vol.trade.benefits.country <- lm(trade_benefits.country ~ volatility, data = data.working)
#Test H1d
vol.trade.support <- lm(trade_support ~ volatility, data = data.working)
vol.new.job <- lm(new_job ~ volatility, data = data.working)
vol.new.job.time <- lm(new_job_time ~ volatility, data = data.working)

stargazer(vol.employment, vol.lose.job, vol.new.job, vol.new.job.time, vol.trade.benefits.self, vol.trade.benefits.region, vol.trade.benefits.country, vol.trade.support, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Main Paper Table 2:
#Test H2a by replicating H1a-H1d with an interaction term with women
#R automatically includes the constitutive terms of the interaction

#Test H1a with women interaction
vol.employment.women <- lm(employment ~ volatility*women, data = data.working)
#Test H1b with women interaction
vol.lose.job.women <- lm(lose_job ~ volatility*women, data = data.working)
#Test H1c with women interaction
vol.trade.benefits.self.women <- lm(trade_benefits.self ~ volatility*women, data = data.working)
vol.trade.benefits.region.women <- lm(trade_benefits.region ~ volatility*women, data = data.working)
vol.trade.benefits.country.women <- lm(trade_benefits.country ~ volatility*women, data = data.working)
#Test H1d with women interaction
vol.trade.support.women <- lm(trade_support ~ volatility*women, data = data.working)
vol.new.job.women <- lm(new_job ~ volatility*women, data = data.working)
vol.new.job.time.women <- lm(new_job_time ~ volatility*women, data = data.working)

stargazer(vol.employment.women, vol.lose.job.women,  vol.new.job.women, vol.new.job.time.women, vol.trade.benefits.self.women, vol.trade.benefits.region.women, vol.trade.benefits.country.women, vol.trade.support.women, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-values for text
summary(vol.new.job.women) # volatility p=0.07
summary(vol.trade.benefits.self.women) # volatility p=0.06
summary(vol.trade.benefits.country.women)# wolatility*women p=0.06

# For descriptive purposes, examine whether women believe they are less likely to find a new job due to trade than men
# Then test how volatility affects the likelihood of finding a new job for working women
women.new.job <- lm(new_job ~ women, data = data.working)
summary(women.new.job)

#Test H2b
vol.new.job.work <- lm(new_job ~ volatility*women, data = data.working)
summary(vol.new.job.work)

# For descriptive purposes, examine whether women believe it will take them longer to find a new job than men
# Then test how volatility affects length of time to find a new job for working women
women.new.job.time <- lm(new_job_time ~ women, data = data.working)
summary(women.new.job.time)
#Test H2c
vol.new.job.time.work <- lm(new_job_time ~ volatility*women, data = data.working)
summary(vol.new.job.time.work)

# To test H2d we first create two data subset for men and women
data.women <- subset(data[data$women==1, ])
data.men <- subset(data[data$gender==1, ])

#Main Paper Table 3
# We then use the subsetted data to repeat H2a-H2c with an interaction with "working"
#Test H1a with "working" interaction on subset of data of women
vol.employment.work <- lm(employment ~ volatility*working, data = data.women)
#Test H1b with women interaction
vol.lose.job.work <- lm(lose_job ~ volatility*working, data = data.women)
#Test H1c with women interaction
vol.trade.benefits.self.work <- lm(trade_benefits.self ~ volatility*working, data = data.women)
vol.trade.benefits.region.work <- lm(trade_benefits.region ~ volatility*working, data = data.women)
vol.trade.benefits.country.work <- lm(trade_benefits.country ~ volatility*working, data = data.women)
#Test H1d with women interaction
vol.trade.support.work <- lm(trade_support ~ volatility*working, data = data.women)
#Test H2b
vol.new.job.work <- lm(new_job ~ volatility*working, data = data.women)
#Test H2c
vol.new.job.time.work <- lm(new_job_time ~ volatility*working, data = data.women)

stargazer(vol.employment.work, vol.lose.job.work, vol.new.job.work, vol.new.job.time.work, vol.trade.benefits.self.work, vol.trade.benefits.region.work, 
          vol.trade.benefits.country.work, vol.trade.support.work, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-value for text
summary(vol.trade.support.work) # volatility p < 0.01



##
## Supplementary Analysis for Appendix
##

# Appendix Table 1 - Study demographics
# Outputs were manually entered into the table in Latex

age1824 <- length(data$age[data$age<=24 & data$age>=18 & is.na(data$age)==FALSE])
age2534 <- length(data$age[data$age<=34 & data$age>=25 & is.na(data$age)==FALSE])
age3544 <- length(data$age[data$age<=44 & data$age>=35 & is.na(data$age)==FALSE])
age4554 <- length(data$age[data$age<=54 & data$age>=45 & is.na(data$age)==FALSE])
age5564 <- length(data$age[data$age<=64 & data$age>=55 & is.na(data$age)==FALSE])
age65up <- length(data$age[data$age>=65 & is.na(data$age)==FALSE])
length(data$age[is.na(data$age)==FALSE]) #5643

age1824/5643 # 8.9%
age2534/5643 # 14.6%
age3544/5643 # 18.6% 
14.6+18.6 #33.2 25-44 year olds
age4554/5643 # 15.8%
age5564/5643 # 20.0%
15.8+ 20 # 35.8 45-64 year olds
age65up/5643 # 22.0%

# Women
table(data$women) 
3000/(3000+2713) #52.5%

# College education
table(data$education)
3363/(3363+2294) #59.4%

# Income
inc0_50 <- length(data$income[data$income==1 | data$income==2 & is.na(data$income)==FALSE]) 
inc51_100 <- length(data$income[data$income==3 | data$income==4 & is.na(data$income)==FALSE])
inc100_200 <- length(data$income[data$income==5 | data$income==6 | data$income==7 | data$income==8 & is.na(data$income)==FALSE])
inc_up <- length(data$income[data$income==9 & is.na(data$income)==FALSE]) 
inc0_50 + inc51_100 + inc100_200 + inc_up #5777

inc0_50/5777 # 37.8%
inc51_100/5777 # 34.2%
inc100_200/5777 #23.8%
inc_up/5777 #4.2%

# Working
table(data$working)
2851/(2851+2862) # 49.9

# white respondents
table(data$white)
4382/(4382+1278)#77.4%

# Import competing - for descriptive statistics
table(data$import_competing)
1227/(1227+4356)

# Appendix Table 2 - Correlation between demographics and treatment assignment
balance <- lm(volatility ~ women + age + education + income + working +white, data = data)
stargazer(balance)

# Appendix Table 3
# Test H3a using interaction term for risk of losing job ("job_risk") 
#Test H1a with job risk interaction
vol.employment.risk <- lm(employment ~ volatility*job_risk, data = data.working)
#Test H1b with job risk interaction
vol.lose.job.risk <- lm(lose_job ~ volatility*job_risk, data = data.working)
#Test H1c with job risk interaction
vol.trade.benefits.self.risk <- lm(trade_benefits.self ~ volatility*job_risk, data = data.working)
vol.trade.benefits.region.risk <- lm(trade_benefits.region ~ volatility*job_risk, data = data.working)
vol.trade.benefits.country.risk <- lm(trade_benefits.country ~ volatility*job_risk, data = data.working)
#Test H1d with job risk interaction
vol.trade.support.risk <- lm(trade_support ~ volatility*job_risk, data = data.working)
#Test H2b
vol.new.job.risk <- lm(new_job ~ volatility*job_risk, data = data.working)
#Test H2c
vol.new.job.time.risk <- lm(new_job_time ~ volatility*job_risk, data = data.working)

stargazer(vol.employment.risk, vol.lose.job.risk, vol.new.job.risk, vol.new.job.time.risk, vol.trade.benefits.self.risk, vol.trade.benefits.region.risk, 
          vol.trade.benefits.country.risk, vol.trade.support.risk, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-values for text
summary(vol.employment.risk) #volatility p=0.03; volatility*job risk p=0.04

# Test whether those doing more unpaid work respond differently to the treatment
# We do this for the full sample, and for the the subset of women
# Using the questions that asked how many hours a week do you spend doing unpaid housework and 
# unpaid care for dependents, we classify those whose answers totalled in the top 1/3 as "high.unpaid.work"
#Test H1a with high.unpaid.work interaction
vol.employment.high.unpaid.work <- lm(employment ~ volatility*high.unpaid.work, data = data)
vol.employment.high.unpaid.work.women <- lm(employment ~ volatility*high.unpaid.work, data = data.women)
#Test H1b with high.unpaid.work interaction
vol.lose.job.high.unpaid.work <- lm(lose_job ~ volatility*high.unpaid.work, data = data)
vol.lose.job.high.unpaid.work.women <- lm(lose_job ~ volatility*high.unpaid.work, data = data.women)
#Test H1c with high.unpaid.work interaction
vol.trade.benefits.self.high.unpaid.work <- lm(trade_benefits.self ~ volatility*high.unpaid.work, data = data)
vol.trade.benefits.self.high.unpaid.work.women <- lm(trade_benefits.self ~ volatility*high.unpaid.work, data = data.women)
vol.trade.benefits.region.high.unpaid.work <- lm(trade_benefits.region ~ volatility*high.unpaid.work, data = data)
vol.trade.benefits.region.high.unpaid.work.women <- lm(trade_benefits.region ~ volatility*high.unpaid.work, data = data.women)
vol.trade.benefits.country.high.unpaid.work <- lm(trade_benefits.country ~ volatility*high.unpaid.work, data = data)
vol.trade.benefits.country.high.unpaid.work.women <- lm(trade_benefits.country ~ volatility*high.unpaid.work, data = data.women)
#Test H1d with high.unpaid.work interaction
vol.trade.support.high.unpaid.work <- lm(trade_support ~ volatility*high.unpaid.work, data = data)
vol.trade.support.high.unpaid.work.women <- lm(trade_support ~ volatility*high.unpaid.work, data = data.women)
#Test H2b
vol.new.job.high.unpaid.work <- lm(new_job ~ volatility*high.unpaid.work, data = data)
vol.new.job.high.unpaid.work.women <- lm(new_job ~ volatility*high.unpaid.work, data = data.women)
#Test H2c
vol.new.job.time.high.unpaid.work <- lm(new_job_time ~ volatility*high.unpaid.work, data = data)
vol.new.job.time.high.unpaid.work.women <- lm(new_job_time ~ volatility*high.unpaid.work, data = data.women)

#Appendix Table 4
stargazer(vol.employment.high.unpaid.work , vol.lose.job.high.unpaid.work ,  vol.new.job.high.unpaid.work , vol.new.job.time.high.unpaid.work , vol.trade.benefits.self.high.unpaid.work , vol.trade.benefits.region.high.unpaid.work , 
          vol.trade.benefits.country.high.unpaid.work , vol.trade.support.high.unpaid.work ,
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 5
stargazer(vol.employment.high.unpaid.work.women , vol.lose.job.high.unpaid.work.women ,  vol.new.job.high.unpaid.work.women , vol.new.job.time.high.unpaid.work.women , vol.trade.benefits.self.high.unpaid.work.women , vol.trade.benefits.region.high.unpaid.work.women , 
          vol.trade.benefits.country.high.unpaid.work.women , vol.trade.support.high.unpaid.work.women , 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))


# Test how provider role based on household income share interacts with volatility
# The measure is coded by quartiles 0-25 being a non-provider, 25-50 ``secondary co-provider'', ``50-75'' ``primary co-provider'' and $>$75\% ``primary or sole provider'' 
# Test interaction term for "provider.role"
#Test H1a with provider role interaction
vol.employment.provider.role <- lm(employment ~ volatility*provider.role, data = data.working)
#Test H1b with provider role interaction
vol.lose.job.provider.role <- lm(lose_job ~ volatility*provider.role, data = data.working)
#Test H1c with provider role interaction
vol.trade.benefits.self.provider.role <- lm(trade_benefits.self ~ volatility*provider.role, data = data.working)
vol.trade.benefits.region.provider.role <- lm(trade_benefits.region ~ volatility*provider.role, data = data.working)
vol.trade.benefits.country.provider.role <- lm(trade_benefits.country ~ volatility*provider.role, data = data.working)
#Test H1d with provider role interaction
vol.trade.support.provider.role <- lm(trade_support ~ volatility*provider.role, data = data.working)
#Test H2b
vol.new.job.provider.role <- lm(new_job ~ volatility*provider.role, data = data.working)
#Test H2c
vol.new.job.time.provider.role <- lm(new_job_time ~ volatility*provider.role, data = data.working)

#Appendix Table 6
stargazer(vol.employment.provider.role , vol.lose.job.provider.role , vol.new.job.provider.role , vol.new.job.time.provider.role , vol.trade.benefits.self.provider.role , vol.trade.benefits.region.provider.role , 
          vol.trade.benefits.country.provider.role , vol.trade.support.provider.role ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 7
# Test H3b using interaction term for willingness to move ("mobility") 
#Test H1a with mobility interaction
vol.employment.mobility <- lm(employment ~ volatility*mobility, data = data.working)
#Test H1b with mobility interaction
vol.lose.job.mobility <- lm(lose_job ~ volatility*mobility, data = data.working)
#Test H1c with mobility interaction
vol.trade.benefits.self.mobility <- lm(trade_benefits.self ~ volatility*mobility, data = data.working)
vol.trade.benefits.region.mobility <- lm(trade_benefits.region ~ volatility*mobility, data = data.working)
vol.trade.benefits.country.mobility <- lm(trade_benefits.country ~ volatility*mobility, data = data.working)
#Test H1d with mobility interaction
vol.trade.support.mobility <- lm(trade_support ~ volatility*mobility, data = data.working)
#Test H2b
vol.new.job.mobility <- lm(new_job ~ volatility*mobility, data = data.working)
#Test H2c
vol.new.job.time.mobility <- lm(new_job_time ~ volatility*mobility, data = data.working)

stargazer(vol.employment.mobility, vol.lose.job.mobility,  vol.new.job.mobility, vol.new.job.time.mobility, vol.trade.benefits.self.mobility, vol.trade.benefits.region.mobility, 
          vol.trade.benefits.country.mobility, vol.trade.support.mobility,
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-value for text
summary(vol.new.job.time.mobility) # volatility*mobility p=.05


#Replicate table 2 for those working Full time (then part time by choice, then part time while searchign for fulltime)

# subset data to only those who are working full time (working full time AND working part time)
data.working.full <- subset(data[data$workforce==1, ])

#Test H1a with women interaction
vol.employment.women.full <- lm(employment ~ volatility*women, data = data.working.full)
#Test H1b with women interaction
vol.lose.job.women.full <- lm(lose_job ~ volatility*women, data = data.working.full)
#Test H1c with women interaction
vol.trade.benefits.self.women.full <- lm(trade_benefits.self ~ volatility*women, data = data.working.full)
vol.trade.benefits.region.women.full <- lm(trade_benefits.region ~ volatility*women, data = data.working.full)
vol.trade.benefits.country.women.full <- lm(trade_benefits.country ~ volatility*women, data = data.working.full)
#Test H1d with women interaction
vol.trade.support.women.full <- lm(trade_support ~ volatility*women, data = data.working.full)
#Test H2b
vol.new.job.women.full <- lm(new_job ~ volatility*women, data = data.working.full)
#Test H2c
vol.new.job.time.women.full <- lm(new_job_time ~ volatility*women, data = data.working.full)

# Using the questions that asked how important a flexible work schedule is, we test whether 
# trade volatility has a different effect for those who place a higher value on a flexible schedule
# We do this for the sample of working respondents, and for the the subset of working women
data.working.women <- subset(data.working[data.working$women==1, ])
#Test H1a with flex.schedule interaction
vol.employment.flex.schedule <- lm(employment ~ volatility*flex.schedule, data = data.working)
vol.employment.flex.schedule.women <- lm(employment ~ volatility*flex.schedule, data = data.working.women)
#Test H1b with flex.schedule interaction
vol.lose.job.flex.schedule <- lm(lose_job ~ volatility*flex.schedule, data = data.working)
vol.lose.job.flex.schedule.women <- lm(lose_job ~ volatility*flex.schedule, data = data.working.women)
#Test H1c with flex.schedule interaction
vol.trade.benefits.self.flex.schedule <- lm(trade_benefits.self ~ volatility*flex.schedule, data = data.working)
vol.trade.benefits.self.flex.schedule.women <- lm(trade_benefits.self ~ volatility*flex.schedule, data = data.working.women)
vol.trade.benefits.region.flex.schedule <- lm(trade_benefits.region ~ volatility*flex.schedule, data = data.working)
vol.trade.benefits.region.flex.schedule.women <- lm(trade_benefits.region ~ volatility*flex.schedule, data = data.working.women)
vol.trade.benefits.country.flex.schedule <- lm(trade_benefits.country ~ volatility*flex.schedule, data = data.working)
vol.trade.benefits.country.flex.schedule.women <- lm(trade_benefits.country ~ volatility*flex.schedule, data = data.working.women)
#Test H1d with flex.schedule interaction
vol.trade.support.flex.schedule <- lm(trade_support ~ volatility*flex.schedule, data = data.working)
vol.trade.support.flex.schedule.women <- lm(trade_support ~ volatility*flex.schedule, data = data.working.women)
#Test H2b
vol.new.job.flex.schedule <- lm(new_job ~ volatility*flex.schedule, data = data.working)
vol.new.job.flex.schedule.women <- lm(new_job ~ volatility*flex.schedule, data = data.working.women)
#Test H2c
vol.new.job.time.flex.schedule <- lm(new_job_time ~ volatility*flex.schedule, data = data.working)
vol.new.job.time.flex.schedule.women <- lm(new_job_time ~ volatility*flex.schedule, data = data.working.women)

#Appendix Table 8
stargazer(vol.employment.flex.schedule , vol.lose.job.flex.schedule , vol.new.job.flex.schedule , vol.new.job.time.flex.schedule , vol.trade.benefits.self.flex.schedule , vol.trade.benefits.region.flex.schedule , 
          vol.trade.benefits.country.flex.schedule , vol.trade.support.flex.schedule ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))
#Appendix Table 9
stargazer(vol.employment.flex.schedule.women , vol.lose.job.flex.schedule.women , vol.new.job.flex.schedule.women , vol.new.job.time.flex.schedule.women , vol.trade.benefits.self.flex.schedule.women , vol.trade.benefits.region.flex.schedule.women , 
          vol.trade.benefits.country.flex.schedule.women , vol.trade.support.flex.schedule.women ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 10
stargazer(vol.employment.women.full , vol.lose.job.women.full , vol.new.job.women.full , vol.new.job.time.women.full , vol.trade.benefits.self.women.full , vol.trade.benefits.region.women.full , 
          vol.trade.benefits.country.women.full , vol.trade.support.women.full ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# subset data to only those who are working part time by choice (working full time AND working part time)
data.working.part.choice <- subset(data[data$workforce==2, ])

#Test H1a with women interaction
vol.employment.women.part.choice <- lm(employment ~ volatility*women, data = data.working.part.choice)
#Test H1b with women interaction
vol.lose.job.women.part.choice <- lm(lose_job ~ volatility*women, data = data.working.part.choice)
#Test H1c with women interaction
vol.trade.benefits.self.women.part.choice <- lm(trade_benefits.self ~ volatility*women, data = data.working.part.choice)
vol.trade.benefits.region.women.part.choice <- lm(trade_benefits.region ~ volatility*women, data = data.working.part.choice)
vol.trade.benefits.country.women.part.choice <- lm(trade_benefits.country ~ volatility*women, data = data.working.part.choice)
#Test H1d with women interaction
vol.trade.support.women.part.choice <- lm(trade_support ~ volatility*women, data = data.working.part.choice)
#Test H2b
vol.new.job.women.part.choice <- lm(new_job ~ volatility*women, data = data.working.part.choice)
#Test H2c
vol.new.job.time.women.part.choice <- lm(new_job_time ~ volatility*women, data = data.working.part.choice)

#Appendix Table 11
stargazer(vol.employment.women.part.choice , vol.lose.job.women.part.choice , vol.new.job.women.part.choice , vol.new.job.time.women.part.choice , vol.trade.benefits.self.women.part.choice , vol.trade.benefits.region.women.part.choice , 
          vol.trade.benefits.country.women.part.choice , vol.trade.support.women.part.choice , 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# subset data to only those who are working part time but looking for full
data.working.part.search <- subset(data[data$workforce==3, ])

#Test H1a with women interaction
vol.employment.women.part.search <- lm(employment ~ volatility*women, data = data.working.part.search)
#Test H1b with women interaction
vol.lose.job.women.part.search <- lm(lose_job ~ volatility*women, data = data.working.part.search)
#Test H1c with women interaction
vol.trade.benefits.self.women.part.search <- lm(trade_benefits.self ~ volatility*women, data = data.working.part.search)
vol.trade.benefits.region.women.part.search <- lm(trade_benefits.region ~ volatility*women, data = data.working.part.search)
vol.trade.benefits.country.women.part.search <- lm(trade_benefits.country ~ volatility*women, data = data.working.part.search)
#Test H1d with women interaction
vol.trade.support.women.part.search <- lm(trade_support ~ volatility*women, data = data.working.part.search)
#Test H2b
vol.new.job.women.part.search <- lm(new_job ~ volatility*women, data = data.working.part.search)
#Test H2c
vol.new.job.time.women.part.search <- lm(new_job_time ~ volatility*women, data = data.working.part.search)

#Appendix Table 12
stargazer(vol.employment.women.part.search , vol.lose.job.women.part.search , vol.trade.benefits.self.women.part.search , vol.trade.benefits.region.women.part.search , 
          vol.trade.benefits.country.women.part.search , vol.trade.support.women.part.search , vol.new.job.women.part.search , vol.new.job.time.women.part.search , 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 13
#Test H1a with "working" interaction on subset of data of men 
vol.employment.work <- lm(employment ~ volatility*working, data = data.men)
#Test H1b with men interaction
vol.lose.job.work <- lm(lose_job ~ volatility*working, data = data.men)
#Test H1c with men interaction
vol.trade.benefits.self.work <- lm(trade_benefits.self ~ volatility*working, data = data.men)
vol.trade.benefits.region.work <- lm(trade_benefits.region ~ volatility*working, data = data.men)
vol.trade.benefits.country.work <- lm(trade_benefits.country ~ volatility*working, data = data.men)
#Test H1d with men interaction
vol.trade.support.work <- lm(trade_support ~ volatility*working, data = data.men)
#Test H2b
vol.new.job.work <- lm(new_job ~ volatility*working, data = data.men)
#Test H2c
vol.new.job.time.work <- lm(new_job_time ~ volatility*working, data = data.men)

stargazer(vol.employment.work, vol.lose.job.work, vol.new.job.work, vol.new.job.time.work, vol.trade.benefits.self.work, vol.trade.benefits.region.work, 
          vol.trade.benefits.country.work, vol.trade.support.work, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# Import competing data subsets 
data.import.comp <- subset(data.working[data.working$import_competing ==1, ]) 
data.NOT.import.comp <- subset(data.working[data.working$import_competing==0, ]) 

#Test H1a with women interaction
vol.employment.import.comp <- lm(employment ~ volatility*women, data = data.import.comp)
vol.employment.NOT.import.com <- lm(employment ~ volatility*women, data = data.NOT.import.comp)
#Test H1b with women interaction
vol.lose.job.import.comp <- lm(lose_job ~ volatility*women, data = data.import.comp)
vol.lose.job.NOT.import.com <- lm(lose_job ~ volatility*women, data = data.NOT.import.comp)
#Test H1c with women interaction
vol.trade.benefits.self.import.comp <- lm(trade_benefits.self ~ volatility*women, data = data.import.comp)
vol.trade.benefits.self.NOT.import.com <- lm(trade_benefits.self ~ volatility*women, data = data.NOT.import.comp)
vol.trade.benefits.region.import.comp <- lm(trade_benefits.region ~ volatility*women, data = data.import.comp)
vol.trade.benefits.region.NOT.import.com <- lm(trade_benefits.region ~ volatility*women, data = data.NOT.import.comp)
vol.trade.benefits.country.import.comp <- lm(trade_benefits.country ~ volatility*women, data = data.import.comp)
vol.trade.benefits.country.NOT.import.com <- lm(trade_benefits.country ~ volatility*women, data = data.NOT.import.comp)
#Test H1d with women interaction
vol.trade.support.import.comp <- lm(trade_support ~ volatility*women, data = data.import.comp)
vol.trade.support.NOT.import.com <- lm(trade_support ~ volatility*women, data = data.NOT.import.comp)
#Test H2b
vol.new.job.import.comp <- lm(new_job ~ volatility*women, data = data.import.comp)
vol.new.job.NOT.import.com <- lm(new_job ~ volatility*women, data = data.NOT.import.comp)
#Test H2c
vol.new.job.time.import.comp <- lm(new_job_time ~ volatility*women, data = data.import.comp)
vol.new.job.time.NOT.import.com <- lm(new_job_time ~ volatility*women, data = data.NOT.import.comp)

#Appendix Table 14
stargazer(vol.employment.import.comp , vol.lose.job.import.comp , vol.new.job.import.comp , vol.new.job.time.import.comp , vol.trade.benefits.self.import.comp , vol.trade.benefits.region.import.comp , 
          vol.trade.benefits.country.import.comp , vol.trade.support.import.comp ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 15
stargazer(vol.employment.NOT.import.com , vol.lose.job.NOT.import.com , vol.new.job.NOT.import.com , vol.new.job.time.NOT.import.com , vol.trade.benefits.self.NOT.import.com , vol.trade.benefits.region.NOT.import.com , 
          vol.trade.benefits.country.NOT.import.com , vol.trade.support.NOT.import.com ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-values for text
summary(vol.trade.support.import.comp) #volatility*women -0.33, p=0.02
# proportion of women supporting trade in import competing for control and treated
t.test(data.import.comp$TRsupp[data.import.comp$women==1 & data.import.comp$volatility==1], data.import.comp$TRsupp[data.import.comp$women==1 & data.import.comp$volatility==0])
#p = 0.03
mean(data.import.comp$TRsupp[data.import.comp$women==1 & data.import.comp$volatility==1 & is.na(data.import.comp$TRsupp)==FALSE])
# 0.49 in treatment
mean(data.import.comp$TRsupp[data.import.comp$women==1 & data.import.comp$volatility==0 & is.na(data.import.comp$TRsupp)==FALSE])
# 0.61 in control


## Effect of volatility based on skill (proxied by education and income) and import competing industry

# Education: We code those as high skilled based on education  
# if they have at least a bachelor's or 4-year college degree
# For income, high income is anyone with household income over $75,000
# We code whether each working respondent is in an import competing industry or not based on census industry codings

#create data subsets based on education, income, and import competing industry
# has a college degree for high education
data.high.edu <- subset(data.working[data.working$education==1, ])
data.low.edu <- subset(data.working[data.working$education==0, ])

# income of $75,000 or greater for high income
data.high.income <- subset(data.working[data.working$income>3, ]) 
data.low.income <- subset(data.working[data.working$income<4, ]) 


# Replicate table 2 for each subset of data
# Education data subsets 
#Test H1a with women interaction
vol.employment.high.edu <- lm(employment ~ volatility*women, data = data.high.edu)
vol.employment.low.edu <- lm(employment ~ volatility*women, data = data.low.edu)
#Test H1b with women interaction
vol.lose.job.high.edu <- lm(lose_job ~ volatility*women, data = data.high.edu)
vol.lose.job.low.edu <- lm(lose_job ~ volatility*women, data = data.low.edu)
#Test H1c with women interaction
vol.trade.benefits.self.high.edu <- lm(trade_benefits.self ~ volatility*women, data = data.high.edu)
vol.trade.benefits.self.low.edu <- lm(trade_benefits.self ~ volatility*women, data = data.low.edu)
vol.trade.benefits.region.high.edu <- lm(trade_benefits.region ~ volatility*women, data = data.high.edu)
vol.trade.benefits.region.low.edu <- lm(trade_benefits.region ~ volatility*women, data = data.low.edu)
vol.trade.benefits.country.high.edu <- lm(trade_benefits.country ~ volatility*women, data = data.high.edu)
vol.trade.benefits.country.low.edu <- lm(trade_benefits.country ~ volatility*women, data = data.low.edu)
#Test H1d with women interaction
vol.trade.support.high.edu <- lm(trade_support ~ volatility*women, data = data.high.edu)
vol.trade.support.low.edu <- lm(trade_support ~ volatility*women, data = data.low.edu)
#Test H2b
vol.new.job.high.edu <- lm(new_job ~ volatility*women, data = data.high.edu)
vol.new.job.low.edu <- lm(new_job ~ volatility*women, data = data.low.edu)
#Test H2c
vol.new.job.time.high.edu <- lm(new_job_time ~ volatility*women, data = data.high.edu)
vol.new.job.time.low.edu <- lm(new_job_time ~ volatility*women, data = data.low.edu)

#Appendix Table 16
stargazer(vol.employment.high.edu , vol.lose.job.high.edu , vol.new.job.high.edu , vol.new.job.time.high.edu , vol.trade.benefits.self.high.edu , vol.trade.benefits.region.high.edu , 
          vol.trade.benefits.country.high.edu , vol.trade.support.high.edu ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 17
stargazer(vol.employment.low.edu , vol.lose.job.low.edu , vol.new.job.low.edu , vol.new.job.time.low.edu , vol.trade.benefits.self.low.edu , vol.trade.benefits.region.low.edu , 
          vol.trade.benefits.country.low.edu , vol.trade.support.low.edu ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# p-values for text
summary(vol.trade.support.low.edu) #volatility*women -0.27 p = 0.04

# Income data subsets 
#Test H1a with women interaction
vol.employment.high.income <- lm(employment ~ volatility*women, data = data.high.income)
vol.employment.low.income <- lm(employment ~ volatility*women, data = data.low.income)
#Test H1b with women interaction
vol.lose.job.high.income <- lm(lose_job ~ volatility*women, data = data.high.income)
vol.lose.job.low.income <- lm(lose_job ~ volatility*women, data = data.low.income)
#Test H1c with women interaction
vol.trade.benefits.self.high.income <- lm(trade_benefits.self ~ volatility*women, data = data.high.income)
vol.trade.benefits.self.low.income <- lm(trade_benefits.self ~ volatility*women, data = data.low.income)
vol.trade.benefits.region.high.income <- lm(trade_benefits.region ~ volatility*women, data = data.high.income)
vol.trade.benefits.region.low.income <- lm(trade_benefits.region ~ volatility*women, data = data.low.income)
vol.trade.benefits.country.high.income <- lm(trade_benefits.country ~ volatility*women, data = data.high.income)
vol.trade.benefits.country.low.income <- lm(trade_benefits.country ~ volatility*women, data = data.low.income)
#Test H1d with women interaction
vol.trade.support.high.income <- lm(trade_support ~ volatility*women, data = data.high.income)
vol.trade.support.low.income <- lm(trade_support ~ volatility*women, data = data.low.income)
#Test H2b
vol.new.job.high.income <- lm(new_job ~ volatility*women, data = data.high.income)
vol.new.job.low.income <- lm(new_job ~ volatility*women, data = data.low.income)
#Test H2c
vol.new.job.time.high.income <- lm(new_job_time ~ volatility*women, data = data.high.income)
vol.new.job.time.low.income <- lm(new_job_time ~ volatility*women, data = data.low.income)
#Appendix Table 18
stargazer(vol.employment.high.income , vol.lose.job.high.income , vol.new.job.high.income , vol.new.job.time.high.income , vol.trade.benefits.self.high.income , vol.trade.benefits.region.high.income , 
          vol.trade.benefits.country.high.income , vol.trade.support.high.income , 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 19
stargazer(vol.employment.low.income , vol.lose.job.low.income , vol.new.job.low.income , vol.new.job.time.low.income ,vol.trade.benefits.self.low.income , vol.trade.benefits.region.low.income , 
          vol.trade.benefits.country.low.income , vol.trade.support.low.income ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))


# Replicate table 2 and for just the US and then the Canadian samples
# First for the US Sample
data.working.US <- subset(data.working[data.working$US_study==1, ])
#Test H1a with women interaction
vol.employment.working.US <- lm(employment ~ volatility*women, data = data.working.US)
#Test H1b with women interaction
vol.lose.job.working.US <- lm(lose_job ~ volatility*women, data = data.working.US)
#Test H1c with women interaction
vol.trade.benefits.self.working.US <- lm(trade_benefits.self ~ volatility*women, data = data.working.US)
vol.trade.benefits.region.working.US <- lm(trade_benefits.region ~ volatility*women, data = data.working.US)
vol.trade.benefits.country.working.US <- lm(trade_benefits.country ~ volatility*women, data = data.working.US)
#Test H1d with women interaction
vol.trade.support.working.US <- lm(trade_support ~ volatility*women, data = data.working.US)
#Test H2b
vol.new.job.working.US <- lm(new_job ~ volatility*women, data = data.working.US)
#Test H2c
vol.new.job.time.working.US <- lm(new_job_time ~ volatility*women, data = data.working.US)

#Appendix Table 20
stargazer(vol.employment.working.US , vol.lose.job.working.US , vol.new.job.working.US , vol.new.job.time.working.US , vol.trade.benefits.self.working.US , vol.trade.benefits.region.working.US , 
          vol.trade.benefits.country.working.US , vol.trade.support.working.US ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# Now for the Canadian sample
data.working.CAN <- subset(data.working[data.working$US_study==0, ])
#Test H1a with women interaction
vol.employment.working.CAN <- lm(employment ~ volatility*women, data = data.working.CAN)
#Test H1b with women interaction
vol.lose.job.working.CAN <- lm(lose_job ~ volatility*women, data = data.working.CAN)
#Test H1c with women interaction
vol.trade.benefits.self.working.CAN <- lm(trade_benefits.self ~ volatility*women, data = data.working.CAN)
vol.trade.benefits.region.working.CAN <- lm(trade_benefits.region ~ volatility*women, data = data.working.CAN)
vol.trade.benefits.country.working.CAN <- lm(trade_benefits.country ~ volatility*women, data = data.working.CAN)
#Test H1d with women interaction
vol.trade.support.working.CAN <- lm(trade_support ~ volatility*women, data = data.working.CAN)
#Test H2b
vol.new.job.working.CAN <- lm(new_job ~ volatility*women, data = data.working.CAN)
#Test H2c
vol.new.job.time.working.CAN <- lm(new_job_time ~ volatility*women, data = data.working.CAN)

#Appendix Table 21
stargazer(vol.employment.working.CAN , vol.lose.job.working.CAN , vol.new.job.working.CAN , vol.new.job.time.working.CAN , vol.trade.benefits.self.working.CAN , vol.trade.benefits.region.working.CAN , 
          vol.trade.benefits.country.working.CAN , vol.trade.support.working.CAN ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# Test whether Covid19's effects on personal employment or industry interact with treatment 
# Indicator variables for respondents who lost their job due to Covid ("covid.job.loss") 
# and those who believe jobs in their industry are in immediate risk due to covid ("covid.job.risk")
#Test H1a with covid job loss and risk interactions
vol.employment.covid.job.loss <- lm(employment ~ volatility*covid.job.loss, data = data.not.working)
vol.employment.covid.job.risk <- lm(employment ~ volatility*covid.job.risk, data = data.working)
#Test H1b with covid job loss and risk interactions
vol.lose.job.covid.job.loss <- lm(lose_job ~ volatility*covid.job.loss, data = data.not.working)
vol.lose.job.covid.job.risk <- lm(lose_job ~ volatility*covid.job.risk, data = data.working)
#Test H1c with covid job loss and risk interactions
vol.trade.benefits.self.covid.job.loss <- lm(trade_benefits.self ~ volatility*covid.job.loss, data = data.not.working)
vol.trade.benefits.region.covid.job.risk <- lm(trade_benefits.region ~ volatility*covid.job.risk, data = data.working)
vol.trade.benefits.country.covid.job.loss <- lm(trade_benefits.country ~ volatility*covid.job.loss, data = data.not.working)
vol.trade.benefits.self.covid.job.risk <- lm(trade_benefits.self ~ volatility*covid.job.risk, data = data.working)
vol.trade.benefits.region.covid.job.loss <- lm(trade_benefits.region ~ volatility*covid.job.loss, data = data.not.working)
vol.trade.benefits.country.covid.job.risk <- lm(trade_benefits.country ~ volatility*covid.job.risk, data = data.working)
#Test H1d with covid job loss and risk interactions
vol.trade.support.covid.job.loss <- lm(trade_support ~ volatility*covid.job.loss, data = data.not.working)
vol.trade.support.covid.job.risk <- lm(trade_support ~ volatility*covid.job.risk, data = data.working)
#Test H2b
vol.new.job.covid.job.loss <- lm(new_job ~ volatility*covid.job.loss, data = data.not.working)
vol.new.job.covid.job.risk <- lm(new_job ~ volatility*covid.job.risk, data = data.working)
#Test H2c
vol.new.job.time.covid.job.loss <- lm(new_job_time ~ volatility*covid.job.loss, data = data.not.working)
vol.new.job.time.covid.job.risk <- lm(new_job_time ~ volatility*covid.job.risk, data = data.working)

#Appendix Table 22
stargazer(vol.employment.covid.job.risk , vol.lose.job.covid.job.risk , vol.new.job.covid.job.risk , vol.new.job.time.covid.job.risk , vol.trade.benefits.self.covid.job.risk , vol.trade.benefits.region.covid.job.risk , 
          vol.trade.benefits.country.covid.job.risk , vol.trade.support.covid.job.risk ,
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Appendix Table 23
stargazer(vol.employment.covid.job.loss , vol.lose.job.covid.job.loss , vol.new.job.covid.job.loss , vol.new.job.time.covid.job.loss, vol.trade.benefits.self.covid.job.loss , vol.trade.benefits.region.covid.job.loss , 
          vol.trade.benefits.country.covid.job.loss , vol.trade.support.covid.job.loss ,  
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

# Replicate table 1 and 2 for the sample of respondents who passed the manipulation check
data.working.M.Check <- subset(data.working[data.working$M.Check==1, ])
#Test H1a
vol.employment.M.Check <- lm(employment ~ volatility, data = data.working.M.Check)
#Test H1b
vol.lose.job.M.Check <- lm(lose_job ~ volatility, data = data.working.M.Check)
#Test H1c
vol.trade.benefits.self.M.Check <- lm(trade_benefits.self ~ volatility, data = data.working.M.Check)
vol.trade.benefits.region.M.Check <- lm(trade_benefits.region ~ volatility, data = data.working.M.Check)
vol.trade.benefits.country.M.Check <- lm(trade_benefits.country ~ volatility, data = data.working.M.Check)
#Test H1d
vol.trade.support.M.Check <- lm(trade_support ~ volatility, data = data.working.M.Check)
vol.new.job.M.Check <- lm(new_job ~ volatility, data = data.working.M.Check)
vol.new.job.time.M.Check <- lm(new_job_time ~ volatility, data = data.working.M.Check)

#Percent of respondents in treatment passing manipulation check
table(data$M.Check[data$volatility==1])
1309/(1309+1525) #46%
#Appendix Table 24
stargazer(vol.employment.M.Check, vol.lose.job.M.Check, vol.new.job.M.Check, vol.new.job.time.M.Check, vol.trade.benefits.self.M.Check, vol.trade.benefits.region.M.Check, vol.trade.benefits.country.M.Check, vol.trade.support.M.Check, 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#Test H2a by replicating H1a-H1d with an interaction term with women
#Test H1a with women interaction
vol.employment.women.M.Check <- lm(employment ~ volatility*women, data = data.working.M.Check)
#Test H1b with women interaction
vol.lose.job.women.M.Check <- lm(lose_job ~ volatility*women, data = data.working.M.Check)
#Test H1c with women interaction
vol.trade.benefits.self.women.M.Check <- lm(trade_benefits.self ~ volatility*women, data = data.working.M.Check)
vol.trade.benefits.region.women.M.Check <- lm(trade_benefits.region ~ volatility*women, data = data.working.M.Check)
vol.trade.benefits.country.women.M.Check <- lm(trade_benefits.country ~ volatility*women, data = data.working.M.Check)
#Test H1d with women interaction
vol.trade.support.women.M.Check <- lm(trade_support ~ volatility*women, data = data.working.M.Check)
#Test H2b
vol.new.job.women.M.Check <- lm(new_job ~ volatility*women, data = data.working.M.Check)
#Test H2c
vol.new.job.time.women.M.Check <- lm(new_job_time ~ volatility*women, data = data.working.M.Check)

#Appendix Table 25
stargazer(vol.employment.women.M.Check , vol.lose.job.women.M.Check , vol.new.job.women.M.Check , vol.new.job.time.women.M.Check ,vol.trade.benefits.self.women.M.Check , vol.trade.benefits.region.women.M.Check , 
          vol.trade.benefits.country.women.M.Check , vol.trade.support.women.M.Check , 
          dep.var.labels = c("Employment", "Lose Job", "New \n Job", "Job Search \n Time", "Self", "Region", "Country", "Trade \n Support"),
          omit.stat = c("rsq", "adj.rsq", "ser", "f"))


# Mediation analysis For Figures 1 - 3 of Appendix
# the gender variable was switched to "women" to be consistent with earlier analysis
#subset for complete observations for mediation analysis
data.working.complete <- data.working[is.na(data.working$trade_support)==FALSE & is.na(data.working$employment)==FALSE &
                                        is.na(data.working$volatility)==FALSE & is.na(data.working$age)==FALSE &
                                        is.na(data.working$education)==FALSE & is.na(data.working$income)==FALSE &
                                        is.na(data.working$women)==FALSE & is.na(data.working$ideology)==FALSE &
                                        is.na(data.working$white)==FALSE & is.na(data.working$import_competing)==FALSE &
                                        is.na(data.working$job_risk)==FALSE & is.na(data.working$mobility)==FALSE &
                                        is.na(data.working$flex.schedule)==FALSE & is.na(data.working$covid.job.loss)==FALSE &
                                        is.na(data.working$covid.job.risk)==FALSE & is.na(data.working$new_job_time)==FALSE &
                                        is.na(data.working$lose_job)==FALSE & is.na(data.working$new_job)==FALSE &
                                        is.na(data.working$wage)==FALSE, ]

# Mediation effect of perceived employment opportunities on trade support
set.seed(1235) # set seed to stabilize p-value
med.employment1 <- lm(employment ~ volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.employment2 <- lm(trade_support ~ employment + volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.employment3 <- mediate(med.employment1, med.employment2, treat="volatility", mediator="employment", sims=1500)
plot(med.employment3) #export size 3.3 x4.3
summary(med.employment3) # ACME p= 0.70

# Mediation effect of estimated time to find a new job on trade support
set.seed(1235) # set seed to stabilize p-value
med.new_job_time1 <- lm(new_job_time ~ volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.new_job_time2 <- lm(trade_support ~ new_job_time + volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.new_job_time3 <- mediate(med.new_job_time1, med.new_job_time2, treat="volatility", mediator="new_job_time", sims=1500)
plot(med.new_job_time3) 
summary(med.new_job_time3) # ACME p=0.064 

## Sensitivity test for mediation analysis on time to find a new job
sens.med <- medsens(med.new_job_time3)
par.orig <- par(mfrow = c(2,2))
plot(sens.med, main="", ylim=c(-1,1))

# Mediation effect of likelihood of losing job on trade support
set.seed(1235) # set seed to stabilize p-value
med.lose_job1 <- lm(lose_job ~ volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.lose_job2 <- lm(trade_support ~ lose_job + volatility + age + education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.lose_job3 <- mediate(med.lose_job1, med.lose_job2, treat="volatility", mediator="lose_job", sims=1500)
plot(med.lose_job3)
summary(med.lose_job3) # ACME p=0.57

# Mediation effect of likelihood of finding new job on trade support
set.seed(1235) # set seed to stabilize p-value
med.new_job1 <- lm(new_job ~ volatility + age +  education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.new_job2 <- lm(trade_support ~ new_job + volatility + age +  education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.new_job3 <- mediate(med.new_job1, med.new_job2, treat="volatility", mediator="new_job", sims=1500)
plot(med.new_job3)
summary(med.new_job3) # ACME p=0.19

# Mediation effect of perceived effect on wages on trade support
set.seed(1235) # set seed to stabilize p-value
med.wage1 <- lm(wage ~ volatility + age +  education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.wage2 <- lm(trade_support ~ wage + volatility + age +  education + income + women + ideology + white + import_competing + job_risk + mobility + flex.schedule + covid.job.loss + covid.job.risk, data=data.working.complete)
med.wage3 <- mediate(med.wage1, med.wage2, treat="volatility", mediator="wage", sims=1500)
plot(med.wage3)
summary(med.wage3) # ACME p=0.36

# Mediation Power analysis
# Mediation power analysis was calculated using https://schoemanna.shinyapps.io/mc_power_med/
# We input the standard deviation for treatment "X" (0.5), the mediator "M" (0.64), and the outcome "Y" (0.93)
# We set the "objective" to vary N to test the power, beghining with an N of 100 and increasing in increments of
# 100 until a maximum N of 4,000.  The simulations were run with 1000 replications and 10,000 Monte Carlo Draws per Rep
# We ran the models for correlations between X, Y, and M of 0.1, 0.075, and 0.05
# Outputs from the simulations were entered into the following CSV file

sim.med.pwr <-read.csv("Mediation_Power_Analysis.csv", header = TRUE)

plot(c(0,4000), c(0,1), type="n", main="Mediation Power Analysis Across Sample Sizes", xlab = "Sample Size", ylab = "Power")
lines(sim.med.pwr$N, sim.med.pwr$Power1, lty = 1, lwd = 2, col = "black")
lines(sim.med.pwr$N, sim.med.pwr$Power2, lty = 2, lwd = 2, col = "black")
lines(sim.med.pwr$N, sim.med.pwr$Power3, lty = 3, lwd = 2, col = "black")
legend(2000,.25, lty=c(1,2,3), lwd=c(1,1,1), col=c("black","black","black"), legend=c("Correlation = 0.1", "Correlation = 0.075", "Correlation = 0.05"))
